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A simple approach for computing unsteady aerodynamic forces from simulated measured strain data is proposed 
in this study. First, the deflection and slope of the structure are computed from the unsteady strain. Velocities and 
accelerations of the structure are computed using the autoregressive moving average model, online parameter 
estimator, low-pass filter, and a least-squares curve fitting method, together with analytical derivatives with respect to 
time. Finally, aerodynamic forces over the wing are computed using modal aerodynamic influence coefficient 
matrices, a rational function approximation, and a time-marching algorithm. A cantilevered rectangular wing is used 
to validate the simple approach. Unsteady aerodynamic forces as well as wing deflections, velocities, accelerations, 
and strains are computed using the CFL3D computational fluid dynamics code and the MSC/NASTRAN finite 
element analysis code; and these CFL3D/NASTRAN-based results are assumed as measured quantities. Computed 
deflections, velocities, accelerations, and unsteady aerodynamic forces are compared with the CFL3D/NASTRAN- 
based results. Computed aerodynamic forces based on lifting-surface theory at subsonic speeds are in good agreement 
with the target aerodynamic forces generated using CFL3D code with the Euler equation. This research demonstrates 
the feasibility of obtaining induced drag and lift forces through the use of distributed sensor technology with measured 


strain data. 


Nomenclature 

[A(s)] = modal aerodynamic influence coefficient matrix in 
Laplace domain 

Ajj = ith mode and jth row coefficient for cosine function 

Cc = chord length at typical section 

[C;] = jth aerodynamic lag term matrix in Roger’s 
approximation at Mach M; 

[Do] = constant term matrix in Roger’s approximation 

[D,] = linear term matrix in Roger’s approximation 

[D>] = quadratic term matrix in Roger’s approximation 

Bij = ith mode and jth row coefficient for sine function 

[E] = state transition matrix 

[G] = damping matrix 

[K] = stiffness matrix 

k = discrete time 

LT = number of aerodynamic lag terms in Roger’s 
approximation 

M = Mach numbers 

[M] = mass matrix 

m = number of reduced frequencies 

{N(s)} =  orthonormalized aerodynamic force vector in 
Laplace domain 

{N}y =  orthonormalized aerodynamic force vector at 
discrete time k 

nm = number of modes 

{Qahk = generalized aerodynamic force vector at discrete 
time k 

dp = dynamic pressure 

{qh = generalized coordinates vector at discrete time k 

{aust = master degree-of-freedom vector at discrete time k 

{que(t)}} = measured master degree-of-freedom vector at 
continuous time f 

{duet} = measured master degree-of-freedom vector at 
discrete time k 

{ds} = _ slave degree-of-freedom vector at discrete time k 
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{due} = static equilibrium position of measured master 

degree-of-freedom vector 
= Laplace variable 

T, = time step for time-marching algorithm 

U = far-field air speed 

{X}, = state vector at discrete time k 

Cay = local maximum unsteady strain 

€vms = root-mean-squared level of strain 

an = strain vector at discrete time k 
ji = ith equivalent viscous damping factor 

Neon = root-mean-squared level of noise 

{y(s)} orthonormalized coordinates vector in Laplace 
domain 

{n} x =  orthonormalized coordinates vector at discrete 
time k 

Kp = reduced frequencies (=(@,C/2U)), where p is 
equal to 1,2,...,m 

0; = ith damping factor 

[®] =  eigenmatrix 

[Dy] =  eigenmatrix corresponds to master degrees of 
freedom 

[Ds] = eigenmatrix corresponds to slave degrees of freedom 

Q; = lag frequencies, where j is equal to 1,2, ..., LT 

Odi = ith aeroelastic damped frequency 

Wnj = ith natural frequency of a structure 

Op = frequencies for computing aerodynamic influence 
coefficient matrices 

[*]" = transpose of a matrix [*] 

[*]7! = inverse of a matrix [*] 

{x} = velocity vector of {x} 

{xt = acceleration vector of {+} 


I. Introduction 


EDUCING fuel consumption for modern aircraft is a goal of the 
National Aeronautics and Space Administration (NASA) 
Aeronautics Research Mission Directorate. This goal can be 
accomplished by reducing airframe weight and aerodynamic drag; 
however, reductions in both for a civil transport aircraft is a challenge 
that may require extensive design changes for optimization and/or 
active controls. In general, the same percentage of weight and drag 
reductions can have a similar effect on fuel savings of a transport 
aircraft [1]. 
Real-time measurement of aerodynamic drag force in flight is an 
essential element for implementing an active drag control technique. 
Two major sources of aerodynamic drag on a business jet and a 
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long-haul transport aircraft at cruise speed are viscous drag and 
induced drag, which are approximately 48-53% (one-half) and 
21-38% (one-third) of the total aerodynamic drag [1]. Moreover, 
induced drag comprises more than 90% of the total drag during 
takeoff for a typical transport aircraft [2]. 

Traditionally, lift forces over the wing are measured using a 
pressure gauge; however, the conventional pressure gauge with its 
associated tubing and cabling can create weight and space limitation 
challenges, and pressure data are available only at the discrete 
location of the gauge. A new method to measure lift forces is needed 
in order to overcome the weight and bulk associated with conventional 
pressure gauges. Development of lightweight distributed sensors is a 
critical technology that can allow continuous monitoring of aerodynamic 
surface shape, dynamic loading, and active control of flexible motion 
and drag. 

Flexible and lightweight optical fibers not only revolutionized 
telecommunications but also altered the sensing world. Optical fibers 
can be used as fiber-optic sensors to measure strain and temperature 
[3]. Fiber-optic sensors have been developed to measure colocated 
strain simultaneously with very high accuracy using fiber Bragg 
gratings (FBGs) [3]. Specifically, the fiber-optic strain sensor (FOSS) 
uses a series of FBGs to obtain measurements at intervals as small as 
every half-inch [4] along a fiber and at frequencies of several 
kilohertz [5]. The ability of FBGs to operate at such high frequencies 
makes them an ideal choice for both static and dynamic aerospace 
applications. The methodology of optically measuring aerodynamic 
forces described by Liu et al. [6] is developed based on beam 
deformation theory. A two-camera videogrammetric system is used 
for optical deformation measurements. The data reduction models for 
extracting the normal force and pitching moment use either the local 
displacement and slope change or the global beam deformation 
profile. 

The availability of wing deflections, velocities, and accelerations 
at all element grid points across the structural finite element (FE) 
model [7,8] will allow engineers to undertake more accurate real-time 
analyses of both internal elastic and inertial forces, as well as external 
aerodynamic forces, at any point on the structure. These force values 
over the entire surface of a structure may also find application in 
structural health monitoring, active flexible motion control, and 
active drag reduction. 

This study focuses on the computation of unsteady aerodynamic 
forces over an entire three-dimensional structure based on measured 
strain information. First, deformations of the entire three-dimensional 
structure are obtained using the two-step approach introduced by Pak 
[7]. Next, velocities and accelerations are computed using an 
autoregressive moving average (ARMA) model, online parameter 
estimator [9], low-pass filter, and a least-squares curve fitting method 
[10], together with analytical derivatives with respect to time. The 
unsteady aerodynamic forces are computed from structural deflections, 
velocities, and accelerations along with linear lifting-surface-based 
modal aerodynamic influence coefficient (AIC) matrices and a rational 
function approximation (RFA). 


II. Mathematical Background 


In this study, external unsteady aerodynamic forces are computed 
from measured strain data. Simulated strain data using CFL3D [ Uy 
MSC/NASTRAN [12] code will be assumed as measured strain data. 
In the first section, deflections and slopes of an entire structure are 
computed from measured strain through the use of the two-step 
approach [7]. Velocities and accelerations of the structure are 
computed in the second section using analytical derivatives with 
respect to time. In the last section, unsteady aerodynamic forces are 
computed in the time domain using the  time-marching 
algorithm [13]. 


A. Computation of Wing Deflection from Measured Strain 


Consider the following structural dynamic governing equations of 
motion as shown in Eq. (1): 


Mig}. + [Giigt. + [Kighe = (Qate dq) 
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where [M], [G], and [K] are mass, damping, and stiffness matrices, 
respectively; and {q}, and {Q,}, are the generalized coordinates and 
aerodynamic force vectors at discrete time k, respectively. 

Out-of-plane deflections along FOSSs can be computed from 
measured unsteady strain data {e}, using a piecewise least-squares 
method, an Akima spline, and a linear assumption, as described in the 
two-step approach [7]. These computed deflections along the fibers 
are combined with an FE model of the structure in order to interpolate 
and extrapolate the deflection and slope of the entire structure 
through the use of the system equivalent reduction and expansion 
process (SEREP) [14]. All of the degrees of freedom (DOFs) in the 
FE model can be rearranged, as shown in Eq. (2): 


—~)am\| _ =) 2 
t= {™} =m =[$" |e 


where {qj}, is the master DOF at discrete time k. In this approach, 
deflections along the FOSS computed from the first step of the two- 
step approach [7] are defined as the master DOF. The remaining 
deflections and slopes over all of the structure are defined as slave 
DOFs at discrete time k, {qs}. In Eq. (2), matrices [By] and [5] are 
eigenmatrices corresponding to master and slave DOFs, respectively; 
and {7}, is the orthonormalized coordinates vector at discrete time k. 
Therefore, Eqs. (3) and (4) are derived from Eq. (2): 


{ause = ult: (3) 


{ste = [Osktax (4) 


In Eq. (3), changing the master DOF at discrete time k {qu}, 
to the corresponding measured value {qy.},, along the FOSS 
gives Eq. (5): 


{duetx = lOultde (5) 


where {qj}; is obtained from the first step of the two-step approach 
[7]. Premultiplying [®y]” to Eq. (5) with matrix inversion gives 
Eq. (6) for computing the orthonormalized coordinates vector at 
discrete time k: 


tac = (Ou) Tu) [ul tamede (6) 


and the generalized coordinates vector {q}, of Eq. (7) is obtained 
from substituting Eq. (6) into Eq. (2): 


{Bi = | Fs ((®y)"[®ul)'[®ul’ faves (7) 


B. Computation of Velocity and Acceleration from Computed 
Wing Deflection 

A simple harmonic motion assumption for the computation of 
wing acceleration works with undamped free vibration problems [8], 
but this assumption cannot handle the heavy damping issues 
associated with aeroelastic oscillation problems. Also, the 
orthonormalized coordinate vector {7}, used for the computation 
of velocities in [8] is not fully decoupled because of coupling between 
structural dynamics and unsteady aerodynamics. 

A new approach for the computations of aeroelastic velocity and 
acceleration is proposed in this study. Velocity and acceleration 
vectors at each sensor location at discrete time k, {qe}, and {que}z, 
of an aeroelastic structural motion are computed using Eq. (8) 
together with analytical derivatives with respect to time: 


{Qme(t)} = {due} + {Sets, cos(@gt) + Bi; snvoae| 


i=1 


(8) 
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where o;(= €;@,;) and @g; are the ith damping factor and damped 
frequency, respectively; and nm is the number of modes. A vector 
{que} represents the static equilibrium position of the unsteady wing 
motion. The coefficients A;; and B;;, j = 1,2, ...,n, for the jth row 
element of the vector can be fitted using a least-squares curve fitting 
technique [9,10]. In this study, o; and @,; are computed using an 
ARMA model, an online parameter estimator, and a sine Butterworth 
low-pass filter [15]. It should be noted in Eq. (8) that o; and w,; are 
estimated; therefore, the least-squares curve fitting in this study is 
based on a linear fitting. From Eq. (8), the velocity and acceleration 
are computed using analytical derivatives with respect to time t. The 
FE model is not used for the computation of {que}. {Quen 
and {Ques k- 

Velocity and acceleration vectors over the entire structure are also 
computed using Eqs. (9) and (10) (SEREP transformation): 


tah =| [leu eud"Meul tine) 


tite =| $e [eu leud"lul ined 0) 


C. Computation of Aerodynamic Force from Wing Deflection, 
Velocity, and Acceleration 

First, modal AIC matrices are computed at Mach number M and 
reduced frequencies x,(=(@,C/2U),p =1,2,...,m) using 
lifting-surface theory: 

[A(k))],  [A(x)], .-- [AG n)] 

where C is the chord length at a typical section, and U is a far-field 
airspeed. These modal AIC matrices can be approximated with 
respect to frequency and Laplace variable s using an RFA. In this 
study, Roger’s approximation [Eq. (11)] is selected for the RFA: 


: LT s[Cj] 
[A(8)] = [Dol + sIDi] + 8*[Da] + 9 
j=l j 


(11) 
Substituting Eq. (2) into Eq. (1) and premultiplying [®]" yields 
Eq. (12): 


[OD] MIP] a}. + [O) 1G}, + [P!TATO]}, 


= [®]"{Oa}, = {Nx (12) 


The orthonormalized aerodynamic force vector {N(s)} in the 
Laplace domain is in Eq. (13): 


EN(s)} = qolA(s)Kin(s)} 
- (‘Daltn + s[D,lin(s)} + 8?[Dalin(s)} 
LX stC ]fm(s)} 
+ a) 


j=l 


(13) 


The time-marching algorithm for the computation of the 
orthonormalized aerodynamic force at discrete time k can be 
summarized as follows [13] in Eqs. (14—21): 


{MN} = dp ((Doltmx + [Pilttie + [Dalte + [Cltxh,) = 14) 


(Be + MB et 


(Xe = LENA 1 + [OLB] 5 


(15) 


where 
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[E] = el4l?« (16) 
Ta 
[o] = [ ella de (17) 
0 
-Q1 0 0 
0 -QI1 0 
[A] = (18) 
0 0 one 
I 
I 
[B]=| . (19) 
I 
[C] = [CiCo, ...,Crz] (20) 
x) 
y=4 (21) 
Xir Jk 


and T,, is a sampling time. Orthonormalized coordinate vectors {};, 
{mp},, and {7}, are computed from Eqs. (22-24): 


tm = (Ou) Tu) [ul tamede (22) 
{pc = (Ou) Tu) [ul tauede (23) 
tie = (Oulu) [ul tauede (24) 


From Eq. (12), the generalized aerodynamic force vector at 
discrete time k, {Qa}, is shown in Eq. (25): 


{Oahe = (QIN, (25) 


FE model 
independent 
: a | 


@e >” 


Measure strain Compute 
along each FOSS aerodynamic load 
4 t 


— (ale ne 


Ope (Bx Ge 
Expand the responses 
to the entire structure 


Compute the slopes and 
deflections along each fiber 


“ Sn #7 FE model 
tach dependent 
(me}x 


Compute accelerations and 
velocities along each fiber 


Fig. 1 Steps used to compute aerodynamic force from measured strain. 
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A rectangular matrix [®]’ can be inverted using a singular-value 
decomposition technique. The steps used to compute aerodynamic 
force from measured strain are depicted in Fig. 1. 

In general, aerodynamic force vectors from a lifting-surface theory 
are normal to the aerodynamic model configuration. Therefore, 
unsteady induced drag force as well as lateral force can be defined 
using the surface normal vector during the unsteady wing surface 
oscillation, as shown in Fig. 2. 


D. Summary of Computation of Aerodynamic Force from Strain 


The following matrices should be calculated before starting the 
computation of unsteady aerodynamic loads during flight: 

1) The first matrices are modal AIC matrices and corresponding 
matrices from RFA: [Do], [D;], [Do], [CyCo, ...,Czr], LE] in 
Eqs. (14) and (16), and [0] [B] in Eq. (15). Matrices [E] and [6] [B] are 
the function of aerodynamic lag frequencies for the Roger’s 
approximation. 

2) The second matrices are the transformation matrices based 
on the SEREP approach and singular value decomposition: 
(®u)' [ul '[®yl’ and ((®]")". 

Step 1) Collect unsteady strain {e},. 
Step 2) Compute wing deflection along the FOSS line, {qy.},, using 
the two-step approach. 


aary Lift from linear panel code 


pour 


a) ZAERO aerodynamic model 


Normal 
ul | 
—> Drag 


Flow 
—> Drag 
Lint Normal 


b) Side view of the unsteady wing motion 


Fig. 2 Definition of the unsteady aerodynamic forces from a linear 
lifting-surface theory. 
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Table 1 Detailed material properties 
of the cantilevered plate wing 

Properties of aluminum insert Number 
Young’s modulus E 9207766 psi 
Shear modulus G 3836570 psi 
Density 0.1 Ib/in.? 
Total weight 0.3806 Ib 
Xcg 2.28 in. 
Ycg 5.75 in. 
Thickness 0.065 in. 


Table 2 Measured and computed natural 


frequencies 
Mode Measured, Hz Computed, Hz Comment 
1 14.29 14.29 First bending 
2 80.41 80.17 First torsion 
3 89.80 89.04 Second bending 


Step 3) Compute wing velocity and acceleration along the FOSS line: 
{duehk and {duet 

Step 4) Compute {7}, {7},, and {7}, using Eqs. (22-24). 

Step 5) Compute the orthonormalized aerodynamic force vector {N}, 
using Eqs. (14) and (15). 

Step 6) Convert the orthonormalized aerodynamic force vector {N}, 
to the generalized aerodynamic force vector {Q,}, using Eq. (25). 
The z-directional load is the lift load. ~ 
Step 7) Compute induced drag and lateral forces using surface normal 
vectors together with the lift force in step 6. 

Steps 1 through 3 are the model independent procedures. On the 
other hand, steps 4 through 7 are dependent on the structural dynamic 
model ® and the unsteady aerodynamic model, [Do], [D;], [D2], 
[C,Co, ..., Cz], [E], and [6] [B]. 


Ill. Results and Discussions 


A cantilevered rectangular wing, shown in Fig. 3, was selected for 
the validation of the proposed approach. This wing, with 6% circular 
arc cross sections and an aspect ratio of 5.0, was built and tested at the 
NASA Langley Research Center (Hampton, Virginia) in 1959 [16]. 
The model had a uniform chord length of 4.56 in., a span length of 
11.5 in., and a thickness of 0.065 in. of aluminum insert covered with 
flexible plastic foam (Fig. 3). The material properties of the aluminum 
insert were assigned a Young’s modulus E of 9.208 Msi 
(Msi = 1,000,000 psi); a shear modulus G of 3.837 Msi; and a 
mass density of 0.1 Ib/in.*. The shaped lumped weights were used to 
match the local cross-sectional weight distribution of the plastic foam. 
Therefore, the small lumped weights were used near the leading and 
trailing edges, and the large lumped weights were used near the 
midchord area. Detailed material properties are shown in Table 1. 
To represent the six FOSSs, the model was fit with 300 beam elements 


SSS = 
SSS SS 
SIE pt 
Ly guaagittrs, — 
|_| E 4 
ee | 
eg KY 
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BEI ] 
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= . ((|| a 
_ il S 
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Flow direction 
Fig. 4 CFD grid for CFL3D computations based on Euler grid. 
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a) M=0.714 
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b) M=0.875 


Fig. 5 Local Mach number distributions at Mach 0.714 and 0.875. 


(50 per each fiber) that computed the axial strain along the length of the 
wing. These six simulated FOSSs were assumed to be zero weight with 
zero stiffness (Fig. 3). 

The frequencies and mode shapes of this cantilevered wing model 
are computed using the MSC/NASTRAN code [12]. Measured and 
computed natural frequencies are summarized in Table 2. Unsteady 
aerodynamic forces as well as wing deflections and velocities are 
computed using the CFL3D code [11]. However, acceleration and 
unsteady strain data are not available from the CFL3D computation; 
therefore, the MSC/NASTRAN code is used to compute target 
acceleration and simulated measured strain data. 

A computational fluid dynamics (CFD) grid configuration for the 
CFL3D computations based on the Euler grid is given in Fig. 4. 
The CFD grid is a multiblock (97 x 73 x 57) grid with H-H topology. 
The time-step size of the unsteady CFL3D computation is 
0.000060515 s, and a total of 10,240 time steps are used in this 
computation. The unsteady aerodynamic theory used in Sec. II.C is 
based on a linear lifting-surface theory: ZAERO code [17]. 
Therefore, a subsonic Mach number should be selected for the 
CFL3D computer simulation to minimize a nonlinear transonic 
effect. Local Mach number distributions under steady-state 
conditions with CFL3D computer simulations are given in Fig. 5. 
In this figure, local Mach number distributions at Mach 0.714 are 
high subsonic conditions. The maximum local Mach number reaches 
the 0.8-0.9 range near the center chord, as shown in Fig. Sa. 
Supersonic subregions are observed in the Mach 0.875 case (that is, 
transonic speed), as shown in Fig. 5b. Therefore, a Mach number of 
0.714 with dynamic pressure of 1.455 psi is selected for the validation 
of the current approach. These CFD-based aerodynamic forces are 
assumed as the target forces in this study. 

The MSC/NASTRAN code was used to calculate unsteady strains 
in this study, and these computed strains are considered as the 
measured strains. For the CFL3D computations, structural mode 
shapes should be provided at the CFD grid points. In this study, the 
structural grid points and the CFD grid points are connected using the 
interpolation elements (“RBE3 element” in MSC/NASTRAN 
terminology) instead of using a surface-splining technique. In the 
CFL3D code, unsteady aerodynamic force vectors are computed at 
the centroids of CFD cells. Therefore, a splining between structural 


Centroid of CFD cell RBE3 elements 
Finite from FEM to CFD 


CFD grids 
FE grids 


RBE3 elements 
from CFD to FEM 


Fig. 6 RBE3 elements between structural grid points and CFD grids 
and centroids (FEM, finite element method). 


grid points and these centroids is also needed for the transient 
response computations with the MSC/NASTRAN code. In this study, 
RBE3 elements are also created between structural grids and these 
centroids of CFD cells, as shown in Fig. 6. It should be noted that the 
well-known numerical problems associated with the Harder and 
Desmarais surface-spline technique [18] can be easily overcome 
through the use of the current technique with RBE3 elements. 

The MSC/NASTRAN modal transient response analysis (solution 
112) with 1024 time steps and a step size of 0.00060515 s is used to 
compute the strains, deflections, velocities, and accelerations. Time 
histories of aerodynamic force vectors at centroids of each CFD cell 
over the upper and lower wing surfaces are computed during 
unsteady CFL3D computation. These unsteady aerodynamic force 
vectors at each time step are converted to the applied force vectors at 
structural grids for the modal transient response analysis. The same 
initial velocity conditions used for the CFL3D computation are also 
used for the modal transient response analysis with the first three 
modes. The structural deflection and velocity values at the leading 
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Fig. 7 Deflection and velocity comparisons using CFL3D and MSC/ 
NASTRAN codes. 
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Fig. 8 Time histories of strain under different levels of random white noise. 
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Fig.9 Summary of the least-squares curve fitting and deflection prediction. 


edge of the wingtip section obtained through the use of the CFL3D 
and MSC/NASTRAN codes are compared in Fig. 7. Excellent 
deflection and velocity matching are observed in this figure. Therefore, 
strain values computed from the MSC/NASTRAN code can be used as 
measured strain values to estimate the unsteady aerodynamic forces 
computed using the CFL3D computer simulation with the Euler 
equation. 

Time histories of strain under different levels of random white 
noise are shown in Fig. 8. Figure 8a shows time histories of strain at 
the leading edge of the wing-root section. Random white noise is 
added to the unsteady strain data to demonstrate the robustness of the 
proposed approach. The strain signal-to-noise ratio (SNR) is defined 
as shown in Eq. (26): 


SNR = 20 X logyy ™ (26) 
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Fig. 10 Time histories of Z deflection under SNR = 0 dB. 
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where €,m; and ym; represent the root-mean-squared (rms) level of 
the unsteady strain and added noise, respectively. In this study, the 
SNRs of 10, 6, and 0 dB are used in the computer simulation. 
Corresponding time histories are shown in Figs. 8b-8d. The local 
SNR (LSNR) is defined as follows: 


Emax 
n 


LSNR = 20 x logo (27) 


rms 


where €max is the local maximum absolute unsteady strain value. In 
Fig. 8d, the LSNR value is the same, with the SNR value near 0.33 s. 
The LSNR value is larger than the SNR value before 0.33 s. The 
LSNR value becomes —9.8 dB near 0.59 s. 

In this study, robustness of the proposed least-squares curve 
fitting method [Eq. (8)] is tested using time histories of unsteady 
strain, shown in Fig. 8. A moving time window of 56 time steps is 
used in this curve fitting, as shown in Fig. 9. The least-squares curve 
fitting method in Eq. (8) is anonlinear fitting problem; however, this 
nonlinear fitting problem becomes a linear problem when the 
damping factors and damped aeroelastic frequencies, o; and w;, are 
provided. In this study, a sine Butterworth low-pass filter [15] witha 
cutoff frequency of 200 Hz is used to estimate reasonable 
frequencies and damping factors from unsteady strain data. The 
number of ARMA coefficients is seven, and the sampling time 
for this online estimator is 0.004236 s (eight steps). In this study, 
a recursive least-squares method based on Bierman’s U-D 
factorization algorithm with a forgetting factor of 0.98 is used as an 
online parameter estimator [9]. Once the fitted coefficients {qy,}, 
Aj;;, and B;; are obtained based on the current 56 time steps, then 
deflections are predicted for the next eight time steps. These eight 
steps correspond to the one sampling period for the online 
parameter estimator. As shown in Fig. 9, the damping factors and 
damped aeroelastic frequencies, o; and w,;, are updated with every 
sampling time step. 

Time histories of Z deflection, velocity, and acceleration under 
0 dB SNR are shown in Figs. 10-12, respectively. The least-squares 
curve fitting starts after the converged damping factors and damped 
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Fig. 13 Time histories of total induced drag force under different levels of random white noise. 
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frequencies are obtained; thus, velocities and accelerations are not 
available until 400 steps (0.2414 s), as shown in Figs. 11 and 12. In 
Figs. 10-12, the solid lines and dashed lines represent target values; 


and corresponding deflection, velocity, and acceleration values 


before (metk {Oueske and {Ouesi) and after Age {Be and {Be) 
using the SEREP transformation, respectively. 
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Fig. 15 Time histories of total lift force under different levels of random white noise. 
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The effect of the piecewise least-squares method for the 
computation of the unsteady deflection can be observed during the 
learning period as shown in Fig. 10. Even with noisy strain data 
(LSNR of 8.7 to 1.6 dB), unsteady deflections are successfully 
obtained. The effect of the SEREP transformation can also be 
observed in Figs. 10—12 as the solid line versus the dashed line. Noise 
in the solid line becomes smoother after the SEREP transformation is 
applied. Finally, the effect of the least-squares curve fitting method in 
Eq. (8) can be clearly observed before and after 0.2414 s, as shown in 
Fig. 10. Noise in unsteady deflection during the learning period is 
drastically reduced after the least-squares curve fitting is used. Wing 
deflection, velocity, and acceleration are in excellent agreement with 
corresponding target values, as shown in Figs. 10-12. The proposed 
least-squares curve fitting method together with the analytical time 
derivatives performs excellently, even with an LSNR of —9.8 dB. 

Modal AIC matrices are computed using the ZAERO code at Mach 
0.714. The ZAERO-based unsteady aerodynamic model configura- 
tion is shown in Fig. 2a. Reduced frequencies of 0.0, 0.006, 0.015, 
0.035, 0.08, 0.13, and 0.26 are selected for this computation. Roger’s 
approximation with four aerodynamic lag terms is used for an RFA of 
these modal AIC matrices. The element-by-element least-squares 
curve fitting with a constraint at the steady-state condition, and a 
reduced frequency of zero, is used in the Roger’s approximation 
procedure. Aerodynamics lag frequencies are 11.81 Hz 
(k = 0.0177), 47.22 Hz (k = 0.0707), 106.2 Hz (ck = 0.1591), and 
188.9 Hz (xk = 0.2829). 

The total induced drag, lateral, and lift forces obtained from the 
current approach under different levels of random white noise are 
compared with the corresponding target aerodynamic forces from 
CFL3D computations in Figs. 13—15. The least-squares curve fitting 
method starts at 0.2414 s in Figs. 13—15. It is interesting that the 
computed forces between times of 0 to 0.2414 s are based on 
unsteady deflection only. Velocities and accelerations are assumed to 
be zero during the learning period, as shown in Figs. 11 and 12. The 
effects of noise can be observed in Fig. 13. Computed total induced 
drag forces with an SNR of 0 dB are the most noisy result, as shown in 
Fig. 13d. 

The wing thickness effects on induced drag and lateral forces 
of 0.0353 and 0.0961 Ibf, respectively, were subtracted from the 
CFD-based target force to have zero force at the steady-state 
condition in Figs. 13 and 14. In general, the current approach based 
on lifting-surface theory gave smaller forces than the target values in 
the cases of lift and lateral forces. The computed induced drag forces 
were in good agreement with the corresponding target drag force, as 
shown in Fig. 13. 

Scaled total induced drag, and lateral and lift forces are shown in 
Fig. 16. In [19], unsteady aerodynamic model tuning of the 
aerostructures test wing 2 for accurate flutter prediction was 
performed at two different flight conditions, and the scaling factors 
obtained were 1.2579 and 1.2719. The average scaling factor of 
1.2649 was multiplied to the aerodynamic forces. When this scaling 
factor was multiplied, the lateral and lift forces were in good 
agreement with the corresponding target values computed using the 
CFL3D code, as shown in Figs. 16b and 16c. Accuracy of the induced 
drag force in Fig. 16a became worse than the previous value. 
Therefore, it could be concluded that the small deviations between 
the current method and the CFL3D code with the Euler equation were 
mainly due to the uncertainty in the lifting-surface aerodynamic 
theory used in this study. 

Recommendations for the practical application of unsteady 
deformation computations were given on page 1071 of [7]. First, it is 
recommended in this study that the lifting-surface theory used for the 
unsteady aerodynamic force computations be upgraded to a method 
based on steady and unsteady CFD computations to improve the 
accuracy of the current proposed approach with transonic aeroelasticity 
or high angle-of-attack flow. Second, an active induced drag control 
system based on the current proposed methodology will be a more 
physics-based approach than the drag control system based on 
measuring fuel flow. Finally, a reduced-order aeroelastic equation of 
motion with a smaller matrix size is recommended for an active 
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flexible motion control system, as well as an active induced drag 
control system. 


IV. Conclusions 


In this study, unsteady aerodynamic forces were computed using 
simulated measured strain data. From unsteady strain information, 
unsteady structural deflections were computed using the two-step 
approach. Unsteady velocities and accelerations were computed using 
an autoregressive moving average model, an online parameter estimator, 
a low-pass filter, and a least-squares curve fitting method, together with 
analytical derivatives with respect to time. The deflections, velocities, 
and accelerations at each sensor location were independent of structural 
and aerodynamic models. The distributed strain data together with the 
current proposed approaches could therefore be used as distributed 
deflection, velocity, and acceleration sensors. 

The general structural deflections, velocities, and accelerations 
were converted to the orthonormalized coordinates to compute 
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orthonormalized aerodynamic force vectors using modal aerodynamic 
influence coefficient matrices. The modal aerodynamic influence 
coefficient matrices were fitted in the Laplace domain using Roger’s 
approximation. Laplace-domain aerodynamics were converted to the 
time domain using a time-marching algorithm. The orthonormalized 
aerodynamic force vectors were transformed to the generalized 
coordinates using pseudomatrix inversion-based on singular-value 
decomposition. Finally, induced drag and lateral forces were obtained 
using surface-normal vectors. In general, computed aerodynamic 
forces based on the lifting-surface theory in subsonic speeds were in 
good agreement with the target aerodynamic forces generated using 
the CFL3D code with the Euler equation. This research demonstrated 
the feasibility of sensing induced drag and lift forces through the use of 
distributed sensor technology, together with the fiber-optic strain 
sensor. Thus, an active induced drag control system could be designed 
using these two computed aerodynamic forces, induced drag and lift, to 
improve the fuel efficiency of an aircraft. 

In this study, interpolation elements (RBE3 elements in MSC/ 
NASTRAN terminology) between structural finite elements grids 
and the computational fluid dynamics grids and centroids were 
successfully incorporated with the unsteady aeroelastic computation 
scheme. The numerical problems often associated with the Harder 
and Desmarais surface-splines technique [18] are thus bypassed 
using the current technique with the RBE3 elements. 

It should be emphasized that the deflection, velocity, and 
acceleration computation based on the proposed least-squares curve 
fitting method are validated with respect to the unsteady strain with a 
LSNR of —9.8 dB. Therefore, the current methodology of computing 
unsteady aerodynamic forces can be applied to the actual flight-test 
data. The most critical technology for the success of the proposed 
approach is the robust online parameter estimator because the least- 
squares curve fitting method depends heavily on aeroelastic system 
frequencies and damping factors. 
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